%% POLICY FUNCTIONs IN LIMITED-COMMITMENT CONTRACT
function  [C, Beta, Eta, K, f, ExRet, mu_v, vol_phi] = limited_optimal_policy(V0, Z, C0, Beta0, Eta0, params, grid_params)
%% parameters
temp = num2cell(params);
[rho, r, ra, alpha, sigma, lambda, M, eta_min, eta_max] = temp{:};

[n_grid, ~, PHI_mat, ~, dphi] = grid_params{:};

PHI = PHI_mat(1,:)';

%% Value Function and derivatives
V = [V0(1);V0;V0(end)];

%Initialize matrices
Vphi_f = zeros(n_grid+2,1);
Vphi_b = zeros(n_grid+2,1);
Vphi_c = zeros(n_grid+2,1);

% First Differences w.r.t. phi
Vphi_f(1:end-1) = (V(2:end) - V(1:end-1))/dphi;
Vphi_b(2:end) = (V(2:end) - V(1:end-1))/dphi;
Vphi_c(2:end-1) = 0.5*(Vphi_f(2:end-1) + Vphi_b(2:end-1));

%Second Differences w.r.t. phi
Vphiphi = (Vphi_f - Vphi_b)/dphi;

% Reduce to original grid
V = V(2:end-1);
Vphi_f = Vphi_f(2:end-1);
Vphi_b = Vphi_b(2:end-1);
Vphi_c = Vphi_c(2:end-1);
Vphiphi = Vphiphi(2:end-1); 

%% Initialize for optimization
C = zeros(n_grid,1);
Beta = zeros(n_grid,1);
Eta = zeros(n_grid, 1);

options = optimoptions('patternsearch',...
            'UseVectorized',true,'Display','off',...
            'ConstraintTolerance', 1e-12,...
            'FunctionTolerance', 1e-12,...
            'MeshTolerance', 1e-12,...
            'MaxIterations', 1000);

%% Find optimizers
parfor iii = 1:n_grid
    x = PHI(iii);
    z = Z(iii);
    
    v = V(iii);
    vx_f = Vphi_f(iii);
    vx_b = Vphi_b(iii);
    vx_c = Vphi_c(iii);
    vxx = Vphiphi(iii);

    c0 = C0(iii);
    beta0 = Beta0(iii);
    eta0 = Eta0(iii);

    Aeq = [];
    Beq = [];

    obj = @(p) limited_fun(p, x, z, v, vx_f, vx_b, vx_c, vxx, params);

    pp = patternsearch(obj,[c0, beta0, eta0],[],[],Aeq,Beq,[0,eta_max*z,eta_max],[M,M,eta_max],[],options)

    C(iii) = pp(1);
    Beta(iii) = pp(2);
    Eta(iii) = pp(3);

end

%% Update controls
ExRet = sigma*Eta.*PHI;
K = ((Beta - Eta.*Z).*C.^rho./(sigma*lambda));

mu_v = ra - (C.^(1-rho))/(1-rho) + 0.5*rho*Beta.^2;

f = C - ExRet.*K;

vol_phi = (PHI.*(1-PHI).*Eta);

end

